load "get_spectrum_normal_significant.ncl"
begin
   expName = (/"bx944","by473"/)
   idir = (/"/g/data/p66/ars599/obs/month/", \
            "/g/data/p66/ars599/work_atm_temp/" /)
   iName = (/"HadISST_sst.1970-2019.n96.nc", \
             "ts_Amon_"+expName(0)+"_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(0)+"a_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(0)+"b_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(0)+"c_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(0)+"d_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(1)+"_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(1)+"a_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(1)+"b_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(1)+"c_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_"+expName(1)+"d_piControl_r1i1p1_0001-0100.nc" /)
   iofile =  addfile(idir(0)+iName(0),"r")
   imx1file =  addfile(idir(1)+iName(1),"r")
   imx2file =  addfile(idir(1)+iName(2),"r")
   imx3file =  addfile(idir(1)+iName(3),"r")
   imx4file =  addfile(idir(1)+iName(4),"r")
   imx5file =  addfile(idir(1)+iName(5),"r")
   imy1file =  addfile(idir(1)+iName(6),"r")
   imy2file =  addfile(idir(1)+iName(7),"r")
   imy3file =  addfile(idir(1)+iName(8),"r")
   imy4file =  addfile(idir(1)+iName(9),"r")
   imy5file =  addfile(idir(1)+iName(10),"r")

; --- indices --- ;
; "nino34" then
   n34locNum = (/-5,5,190,240/) ; sLat nLat wLon eLon 
; "atn" then
   atnlocNum = (/-3,3,340,360/) ; sLat nLat wLon eLon
; "nta" then
   ntalocNum = (/-3,3,300,345/) ; sLat nLat wLon eLon
;
; a) pacemaker n34 + obs
; b) pacemaker atn + obs
; c) control atn + obs
;
   obs_sst = iofile->sst
   modx1_sst = imx1file->ts
   modx2_sst = imx2file->ts
   modx3_sst = imx3file->ts
   modx4_sst = imx4file->ts
   modx5_sst = imx5file->ts
   mody1_sst = imy1file->ts
   mody2_sst = imy2file->ts
   mody3_sst = imy3file->ts
   mody4_sst = imy4file->ts
   mody5_sst = imy5file->ts
; n34 obs + pace
   obs_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( obs_sst(:,{ n34locNum(0):n34locNum(1) },{ n34locNum(2):n34locNum(3) }) ),(/1,2/))
   mody1_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mody1_sst(:,{ n34locNum(0):n34locNum(1) },{ n34locNum(2):n34locNum(3) }) ),(/1,2/))
   mody2_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mody2_sst(:,{ n34locNum(0):n34locNum(1) },{ n34locNum(2):n34locNum(3) }) ),(/1,2/))
   mody3_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mody3_sst(:,{ n34locNum(0):n34locNum(1) },{ n34locNum(2):n34locNum(3) }) ),(/1,2/))
   mody4_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mody4_sst(:,{ n34locNum(0):n34locNum(1) },{ n34locNum(2):n34locNum(3) }) ),(/1,2/))
   mody5_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mody5_sst(:,{ n34locNum(0):n34locNum(1) },{ n34locNum(2):n34locNum(3) }) ),(/1,2/))
; atn obs + pace
   obs_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( obs_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   mody1_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( mody1_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   mody2_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( mody2_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   mody3_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( mody3_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   mody4_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( mody4_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   mody5_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( mody5_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
; atn obs + control
   modx1_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( modx1_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   modx2_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( modx2_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   modx3_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( modx3_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   modx4_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( modx4_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))
   modx5_atn = dim_avg_n_Wrap(rmMonAnnCycTLL( modx5_sst(:,{ atnlocNum(0):atnlocNum(1) },{ atnlocNum(2):atnlocNum(3) }) ),(/1,2/))

   obs_n34 = dtrend(obs_n34,False)
   obs_std = dim_stddev( obs_n34 )
   obs_n34 = obs_n34 / obs_std
   obs_n34@long_name = ""

   mody1_n34 = dtrend(mody1_n34,False)
   mody2_n34 = dtrend(mody2_n34,False)
   mody3_n34 = dtrend(mody3_n34,False)
   mody4_n34 = dtrend(mody4_n34,False)
   mody5_n34 = dtrend(mody5_n34,False)

   modx1_atn = dtrend(modx1_atn,False)
   modx2_atn = dtrend(modx2_atn,False)
   modx3_atn = dtrend(modx3_atn,False)
   modx4_atn = dtrend(modx4_atn,False)
   modx5_atn = dtrend(modx5_atn,False)

   mody1_atn = dtrend(mody1_atn,False)
   mody2_atn = dtrend(mody2_atn,False)
   mody3_atn = dtrend(mody3_atn,False)
   mody4_atn = dtrend(mody4_atn,False)
   mody5_atn = dtrend(mody5_atn,False)

   mody_n34 = mody1_n34
   mody_n34 = (mody1_n34+mody2_n34+mody3_n34+mody4_n34+mody5_n34)/5
   mody_std = dim_stddev( mody_n34 )
   mody_n34 = mody_n34 / mody_std
   mody_n34@long_name = ""

   obs_atn = dtrend(obs_atn,False)
   obs_std = dim_stddev( obs_atn )
   obs_atn = obs_atn / obs_std
   obs_atn@long_name = "Observation ATN index"

   mody_atn = mody1_atn
   mody_atn = (mody1_atn+mody2_atn+mody3_atn+mody4_atn+mody5_atn)/5
   mody_std = dim_stddev( mody_atn )
   mody_atn = mody_atn / mody_std
   mody_atn@long_name = "Pacemaker run mean ATN index"

   modx_atn = modx1_atn
   modx_atn = (modx1_atn+modx2_atn+modx3_atn+modx4_atn+modx5_atn)
   modx_std = dim_stddev( modx_atn )
   modx_atn = modx_atn / modx_std
   modx_atn@long_name = "Control run mean ATN index"

;--- output variables ---
  diro    = "./"
  pltName = "fig04"
  outFile = diro+pltName+".nc"
  system("\rm "+outFile)
  fout         = addfile( outFile,"c")
  fout->obs_atn = obs_atn
  fout->modx_atn = modx_atn
  fout->mody_atn = mody_atn


;--- time axial ---;
  time  = mody_n34&time
  tt    = cd_calendar(time, 4)
  wks_time = tt

;*** ploting spectrum ********************
  optPlt          = True
  optPlt@trYMinF = 0.
  optPlt@trXMinF = 0.1
  optPlt@vpHeightF= 0.1          ; change aspect ratio of plot
  optPlt@vpWidthF = 0.9
  optPlt@prdGap   = 1./12
  optPlt@trXLog   = False
  optPlt@frqOrPrd = False
  optPlt@vpHeightF= 0.4                    ; change aspect ratio of plot
  optPlt@vpWidthF = 0.8
;  optPlt@tiMainString = "Power Spectrum of Nino3.4"
  optPlt@units = "Years"
  pltType  = "png"
  pltName  = "fig04_allindices_spec"

    wks  = gsn_open_wks(pltType,pltName)       ; specifies a ps plot
    plot = new(4,graphic)
    optPlt@trYMaxF = 40.
    optPlt@trXMaxF = 12.
  optPlt@gsnCenterString = "Nino3.4 power spectrum of pacemakers"
    plot(0) = plot_spectrum_normal(wks,mody_n34,optPlt)              ; raw data
  optPlt@gsnCenterString = "ATN power spectrum of pacemakers"
    plot(1) = plot_spectrum_normal(wks,mody_atn,optPlt)              ; raw data
  optPlt@gsnCenterString = "ATN power spectrum of base runs"
    plot(2) = plot_spectrum_normal(wks,modx_atn,optPlt)              ; raw data
  optPlt@gsnCenterString = "ATN power spectrum of observation"
    plot(3) = plot_spectrum_normal(wks,obs_atn,optPlt)              ; raw data
; --- add legend --- ;
  gres = True
  gres@YPosPercent = 80.    ; expressed as %, 0->100, sets position of top border of legend 
                            ;  when gres@Position is set to its default setting of "Top" (Default = 95.)
  gres@XPosPercent = 60     ; expressed as %, 0->100, sets position of left border of legend(Default = 5.)
  lineres = True
  lineres@lgLineColors = (/"black", "red","orange","dodgerblue"/) ; line colors
  lineres@lgLineThicknesses = 4                          ; line thicknesses
  lineres@LineLengthPercent = 9.                         ; expressed as %, 0->100, length of line
  textres = True
  textres@lgLabels = (/"Pacemaker run"/)  ; legend labels (required)
  plot(0) = simple_legend(wks,plot(0),gres,lineres,textres)
  textres@lgLabels = (/"Pacemaker run"/)  ; legend labels (required)
  plot(1) = simple_legend(wks,plot(1),gres,lineres,textres)
  textres@lgLabels = (/"Control run"/)  ; legend labels (required)
  plot(2) = simple_legend(wks,plot(2),gres,lineres,textres)
  textres@lgLabels = (/"HadISST"/)  ; legend labels (required)
  plot(3) = simple_legend(wks,plot(3),gres,lineres,textres)
; --- panel plot --- ;
    resP = True
;    resP@gsnPanelMainString        = "Power Spectrum of WWB" ; new resource added in NCL V6.4.0
    resP@gsnPanelFigureStringsJust = "TopLeft"                      ; new resource added in NCL V6.4.0
    resP@gsnPanelFigureStrings     = (/"a)","b)","c)","d)"/) ; add strings to panel
    resP@gsnPanelRowSpec = False                   ; tell panel whatorder to plot
    gsn_panel(wks,plot,(/3,1/),resP)
delete(wks)   ; Make sure PS file is closed

end
